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Abstract 

The discrepancy between structural and functional connectivity in neural systems forms the challenge in understanding 
general brain functioning. To pinpoint a mapping between structure and function, we investigated the effects of 
(in)homogeneity in coupling structure and delays on synchronization behavior in networks of oscillatory neural masses by 
deriving the phase dynamics of these generic networks. For homogeneous delays, the structural coupling matrix is largely 
preserved in the coupling between phases, resulting in clustered stationary phase distributions. Accordingly, we found only 
a small number of synchronized groups in the network. Distributed delays, by contrast, introduce inhomogeneity in the 
phase coupling so that clustered stationary phase distributions no longer exist. The effect of distributed delays mimicked 
that of structural inhomogeneity. Hence, we argue that phase (de-)synchronization patterns caused by inhomogeneous 
coupling cannot be distinguished from those caused by distributed delays, at least not by the naked eye. The here-derived 
analytical expression for the effective coupling between phases as a function of structural coupling constitutes a direct 
relationship between structural and functional connectivity. Structural connectivity constrains synchronizability that may be 
modified by the delay distribution. This explains why structural and functional connectivity bear much resemblance albeit 
not a one-to-one correspondence. We illustrate this in the context of resting-state activity, using the anatomical 
connectivity structure reported by Hagmann and others. 
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Introduction 

Much of the current focus in the empirical study of large-scale 
neuronal networks has been on their intrinsic activity and the 
degree to which the coherent patterns of this intrinsic activity 
reflect anatomy. The use of fMRI and diffusion spectrum imaging 
has allowed for a comprehensive evaluation of the structure- 
function map of resting-state networks (RSNs). In fMRI the spatial 
patterns of spontaneous changes in blood oxygenation level- 
dependent signals seem to reflect the generating neural architec- 
ture of RSNs. Despite the very slow changes of these signals, 
Biswal and co-workers [1] defined RSNs as networks of brain 
areas that exhibit temporally coherent activity in the absence of 
identifiable externally imposed or measurable events. More 
recently, RSNs penetrated the field of encephalography [2,3]. 
For M/EEG, locally synchronized neural activity is considered to 
yield macroscopic oscillations that provide a basis for defining 
functional brain networks [4]. In most studies, structural 
connectivity is considered a good predictor of functional connec- 
tivity [5,6]: Structural connectivity agrees with the anatomical 
connections between network nodes and functional connectivity 
covers the statistical relationship of nodal activity. 

The predictive value of structure for function found support in 
recent modeling work using full brain systems with realistic 
anatomy, which demonstrated the structural dependency of 



functional network configurations [7]. There, functional connec- 
tivity has been estimated between all nodes over several hundred 
seconds of simulated time yielding the pattern of functional 
connectivity over this time window that largely reproduced the 
structural connectivity. At smaller time windows, however, 
shorter-living patterns of functional connectivity emerged that 
had not been predicted by anatomy. To understand this 
discrepancy we investigated effects of time delays vis-a-vis effects 
of structural inhomogeneity on synchronization patterns of 
neuronal networks. 

Delays are inherent in neuronal networks due to finite 
conduction velocities [8] and synaptic transmission [9]. Ignoring 
delays may be a valid starting point for mathematical analysis but 
when doing so one runs the risk of loosing biological plausibility. 
However, incorporating delays in oscillatory networks does come 
with immense challenges. Already for low-dimensional oscillatory 
systems (or for high-dimensional ones with strong symmetry) the 
presence of delays is known to change the dynamical repertoire 
significandy [10,11]. Yeung and Strogatz showed for very large 
networks how time delays can alter synchronization properties, 
even if the structure is isotropic and homogeneous [12]; see also 
[13,14]. Numerical assessments revealed similar results for 
biologically motivated and hence more inhomogeneous connec- 
tivities. Delays seem to be crucial in establishing the spatio- 
temporally organized fluctuations typically observed in resting 
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Author Summary 

Separating the time scale of oscillations from that of the 
phase dynamics allowed for reducing a network of 
coupled neural mass models to a system of phase 
oscillators. We studied the dynamics of networks of phases 
and their synchronization characteristics as being seminal 
for functional neural networks. We put particular focus on 
effects of time delays in the coupling on the network 
dynamics and contrasted that to effects due to altered 
structural connectivity. Does neuroanatomical structure 
prescribe all the macroscopic activity patterns that we 
observe through electrophysiological brain recordings? We 
found that heterogeneity in structural coupling and 
distributed delays have equivalent effects on the shape 
of phase distributions, i.e., on functional connectivity. The 
contribution of changes in structural connectivity to 
network synchronization can therefore not readily be 
distinguished from that of distributed delays. Interestingly, 
the emergence of phase clusters in networks requires a 
subtle interplay between coupling and delays, which may 
form a window into disentangling structural effects from 
those induced by delay distributions. Therefore, when 
investigating neural network behavior, both structural 
connectivity and delay distribution should be addressed. 



state brain recordings [15-17]. In the present study we sought to 
tackle this issue and separated the effect of time delays from that of 
inhomogeneous connectivity by studying networks consisting of 
distinct neural masses. Neural mass models offer a low-dimensional 
description of the dynamics of a large neuronal population and exist 
in a variety of forms [18]. We chose for Freeman's seminal model 
[19-21], since it covers the dynamics of mean membrane potential 
changes that relate closely to encephalographic signals. A network of 
such entities may constitute RSNs if we regard the neural masses to 
be representative of individual brain areas. 

Throughout the paper we describe functional connectivity by 
means of phase synchronization whose dynamics can be estimated 
in voltage-based and firing-rate models using a combination of 
rotating wave and slowly varying amplitude approximations, or in 
brief averaging, see [22,23]. In the Methods section this combina- 
tion of approximations is briefly summarized for Freeman neural 
mass models in the oscillatory regime. Central outcome measure is 
thus the phase dynamics of the individual nodes in the network or, 
to be more precise, the density of the nodes' phases as a function of 
time, often also referred as time-dependent population distributions. 
We note that we applied this approach before to instantaneously 
coupled Wilson-Cowan firing rate models [24] (see also [25]) but, as 
said, we here chose for the Freeman model for an easier comparison 
with M/EEG studies. For coupled Freeman models we could 
analytically determine the corresponding stationary distributions 
even in the presence of delays and inhomogeneous coupling 
between neural masses. We could not only prove the existence of 
these solutions, but we were also able to determine the loss of 
stability of the desynchronized state as soon as the overall coupling 
strength exceeded a critical value. More complicated scenarios 
including biological plausible anatomical adjacencies were treated 
numerically to illustrate the non-trivial relationship between 
structural and functional connectivity. 

Results 

We considered a set of N coupled neural masses whose mean 
membrane potentials V k follow the dynamics 



V k = - (at + P k ) V k - a k ji k V k + a k $ k q k + 

N 

<t k p k Y, c ki[w{t)*S[V,{t-x kl )\\ . (1) 

In this expression a k and p k represent mean rise and decay times of 
neural responses in population k, q k stands for an external input, and 
S[-] denotes a sigmoidal activity function covering the effects of pulse- 
coupled neurons in populations 1=1 ...N [19,26]. Corresponding 
mean activities V; arrive at population k after yet arbitrary delays z k i. 
The structural connectivity matrix C k i served to introduce both 
excitatory and inhibitory connections in the thus asymmeuic coupling; 
see Fig 1A. We first considered the case in which a large degree of 
homogeneity was present in C k i to define a 'baseline'. Subsequently we 
introduced inhomogeneity to mimic, e.g., the sparse connectivity 
presumably underlying RSNs. Two seminal coupling schemes are 
sketched in Fig 1 . Excitatory and inhibitory populations were always 
properly balanced to stabilize oscillatory behavior [27,28]. This 
translates to the condition that at least one pair of the eigenvalues of 
the linearized system around the fixed points Vjp,Uj^ must be 
imaginary with positive real part. From the Methods section it follows 
that the considered coupling schemes did satisfy this condition. 

For the convolution we chose an exponentially decaying 
kernel W 



W{t) 



ye~ yt for t>0 
0 otherwise ' 



(2) 



to represent the dynamics at the synaptic junction. In the particular 
case of infinitesimal memory, i.e. for y— >co, we found the 
corresponding phase dynamics by transforming the neural mass 
dynamics to polar coordinates around an unstable focus, i.e. 
V k -» + R k cos(fi? + (j) k ); R k denotes the amplitude of oscillation 
and <j) k its phase corresponding to the central frequency Q. 
Subsequently, we averaged the dynamics over one period 2n/Q. 
We assumed that the characteristic times of amplitude and phase 




Figure 1. Structural connectivity in the case of two coupled 
pairs of excitatory/inhibitory neural masses. A: Homogeneous 
coupling with coupling matrix c* om >. " : Inhomogenous coupling using 
the cou pling matrix C* ,h) . In both cases excitatory/inhibitory-pairs have 
a symmetric 'internal' coupling with strength ±Ki but the coupling 
between these pairs differs. In A the between-pair coupling is 
homogeneous at strength Ki whilst in B the between-pair coupling 
may differ across pairs; i.e. it is randomly chosen from a certain 
distribution U with the constraint that inhibitory units map (on average) 
with negative and excitatory with positive coupling strength. See text 
for more details and Figs 7 and 8A for the more complicated coupling 
schemes employed below. 
doi:10.1371/journal.pcbi.1003736.g001 
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dynamics are significantly larger than the period 2n / SI, i.e. amplitude 
and phase dynamics are slow compared to the oscillation. We also 
assumed the time delays to be of the same order of magnitude or 
smaller than the period. As a result the delays Zkl reduced to mere 
phase shifts A# and the phase dynamics became 



N 

k = (Ok- a k p k 2J D ki sin(<^ - <j> k + A«) 
1=1 

1 Ri 

7 



with 



D kr - 



2D.R k 



•sy?]c kl 



(3) 



In the dynamics (3) the frequencies obeyed the form 



and the phase shifts read A/ f / = 



-Q-iki. As said, 



20. «— • ~- 2 

delays could be transformed to phase shifts due to a time scale 
separation in the system, i.e. the phase dynamics was set sufficiently 
slow compared to the oscillation and the delays %kj were up to the 
same order of magnitude as one period of oscillation (cf. Fig 2). 

At first glance the dynamics (3) seemed to largely resemble the 
Kuramoto network of phase oscillators that is well known for its 
bifurcation scheme from desynchronized to synchronized states. 
The latter, i.e. the fully synchronized state, would imply large-scale 
— if not whole-brain — synchronization, which, apart for very 
pathological cases, is not observed experimentally. Nonetheless we 
considered linking our study to Kuramoto's profound work to be 
very insightful, as understanding the Kuramoto network is 
essential for understanding synchronizability in our more general 
framework. A closer look revealed that, although similar, (3) differs 
from the Kuramoto network in two important aspects. First, in the 
phase dynamics the coupling is given by Dki, which in general is 
not entirely homogeneous as in Kuramoto's case. This expression 
for Dki agrees with the previously derived phase dynamics of 
coupled Wilson-Cowan oscillators in that the amplitude relation 
R/ 1 R k between nodes / and k affects the corresponding (relative) 
phase [24]. Second, the finite delay z k i yields a non-trivial phase 

7T 

shift Aki= — —SlZkl- This phase shift might alter the spectrum of 

phase synchronization entirely; see, e.g., [29,30] for the case. 

Instead of studying the mere collection of phases {^i, . . . ,(j) N } 
we investigated the dynamics of their probability density 



V = V(<j>\, ■ . ■ ,<l>N,t), because it forms a proper measure of 
synchronization. To explain: Synchronization around a certain 
phase value 9 manifests itself as a peak in the probability density 
around that value, i.e. a phase cluster around 6. For our set of 
Freeman neural masses we found that the phase distribution 
follows the dynamics 



w k — ct. k 



Pk D ki sin (fa ~ fa + At/) 



V. (4) 



"Pstat = Tj Y^k-l ^fa ~ ® k ^' because the network has a countable 



This is the continuity equation, which equals the zero-noise limit 
of the corresponding Fokker-Planck equation when presuming 
ergodicity and very large networks (N— > oo). Despite the presence of 
delays ikl we here succeeded to specify stationary solutions 
"Pstat (<^i> • • • ^n)- Crucial in finding these solutions 
was the fact that they can always be written as a mere sum of 
distinct clustered states, i.e. they always obey the form 
1 

~N 2-^i k =\ 

number of nodes. Put differently, the phase distribution contains 
clusters centered around 9i,02, - ■ • In general, the real number 
of clusters is given by the number M of distinguishable centroids 9 k; 
we always consider M<N. The number M strongly depends on 
the distribution of delays and/ or inhomogeneity of the coupling 
matrix Cu or strictly speaking of the effective coupling D ki- 
ln order to illustrate stationary solutions of (4) we first specified the 
coupling matrix Cki to be either homogeneous or inhomogeneous. In 
both cases we labeled excitatory and inhibitory populations by £ and 
I, respectively, and defined the corresponding sets 

AO . „ f , N 

• + !,. 



£=lk=\, 



N > . As mentioned 



, , and X- ... 
2 J I 2 

above we guaranteed throughout analysis that all excitatory 

populations had an inhibitory counterpart generating and stabilizing 

oscillatory activity at around the central frequency O. [27,28]. For the 

sake of simplicity we ignored self-coupling of neuronal populations. 

The stationary phase distribution could thus be written as 



Pstat=P stllt (^A) = ^^(^-0A-)+^E' 5 (^-^) • ( 5 ) 

Y keS k<El 




Figure 2. Frequency and period as a function of delay. A: Oscillation frequency Q. B: Period 2n/Q. The analytical estimate under small-delay 
approximation (see Methods) is displayed as a solid line; results of numerical simulations of the dynamics (1) as dots ( >). Simulations were performed 
for homogeneous coupling C*° m> with within-pair coupling K\ = 1, between-pair coupling strength ^2 = 0.2, and distributed delays t*j; see text for 
details. The presence of between-pair coupling K2 and even heterogeneity in coupling, CjJ om) ->c5' / nh, : did not yield qualitatively different results. 
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Figure 3. Phase distributions for strong coupling (K 2 =0.6) and equal delays. A: Vanishing delay (t = 0). B: Fixed finite delay r-0.05. C: 
Fixed finite delay t = 0.2. In all cases synchronized solutions V stsl t,2 can be seen and appeared to be the only stable solution in the case of strong 
coupling ^2=0.6. Phase distributions were obtained from simulations of the original system (1), where phases were extracted using the Hilbert 
transforms of Vt(t); see Methods section. In this and all subsequent figures excitatory phases (j> £ are depicted in orange and inhibitory phases <j> x in 
green. Radial axes are normalized. 
doi:1 0.1 371 /journal.pcbi.1 003736.g003 



Homogeneous coupling 

In the homogeneous case we employed the coupling scheme 
sketched in Fig 1A. Excitatory and inhibitory populations were fully 
comiected (apart from self-coupling) with coupling values K\ and K 2 
discriminating within-pair and between-pair coupling. For the 
numerical assessment we always fixed the within-pair coupling to 
= 1 . In more detail, we chose the overall homogeneous coupling 
matrix as 



r 



(hom) 



i(hom) 
ff 



c 



(horn) 



SI 



c 



r(hom) 



C 



(horn) 



with sub-population connectivities 



.--i(hom) 
L ff 



and 



2K 2 

\-lsT 



c 



((horn) 



If 



2K 2 

If 



A', 

2K 2 
L ~N~ 



A'i 



.(horn) 



c 



i(hom) 

ex 



In the absence of delays, i.e. for Tjy = 0 (Ay = — ) and sufBciendy strong 

coupling K 2 we found robust distributions with M = 2 phase clusters, 
one containing all the excitatory populations and one all the inhibitory 




Figure 4. Phase distributions for weak coupling (K2 = 0.2) and equal delays. A: Vanishing delay (t = 0). B: Fixed finite delay t= 0.05. C: Fixed 
finite delay t = 0.2. As for strong coupling we found synchronized solutions V st!it 2 in the weak coupling case K2 = 0.2. In addition solutions V s mA 
were also present. The number of clusters did depend on the delay value x. That is, altering the delay from r = 0 to x = 0.2 caused a switch in stability 
between these two solutions. Since for t = 0.05 the two clusters within an £/X-group only showed a small phase difference, we conjecture that 
t = 0.05 is close to the critical value of this bifurcation parameter. 
doi:1 0.1 371 /journal.pcbi.1 003736.g004 
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Figure 5. Phase distributions for different coupling strengths and distributed delays. A: A 2 = 0.6, xu ~W(0,0.05). B: A 2 = 0.6, 
Hz ~t/(0,0.2). C: i^2 = 0.2, n/~W(0,0.05). The effect of a delay distribution and consequently the presence of M=N centroid phase values Ot 
manifested itself as a widening of the phase distribution compared to the constant delay cases in Fig 4. Narrow delay distributions with 
tkl ~ZV(0,0.05), i.e. the th were randomly drawn from a uniform distribution over the interval [0,0.05], yielded comparably narrow phase distributions 
located around two narrow peaks (peaks representing the £- and X-populations; (A)). Increasing the width of the delay distribution to n/~W(0,0.2) 
(B) had a very similar effect as lowering the coupling strength A2 = 0.2 (Q: in both cases the phase distribution widened substantially. 
doi:10.1371/journal.pcbi.1003736.g005 



An example of this solution is illustrated in Fig 3A; the remaining 
panels in that figure refer to cases of non-vanishing delay that will be 
summarized below. We note that due to symmetry the homogeneous 
case with = 0 can be readily transformed, proving its resemblance 
with the Kuramoto network [31]. The somewhat lengthy analytic 
derivations are given in the Methods section. 

Next to the homogeneously synchronized state "P s t a t,2 we found 
a solution with M = 4 phase clusters given by 

PsvMhA) = J - Osj.) + 4 <5(<^ - Qsz) + 



Again we refer to the Methods for the analytical treatment with 
respect to the existence of this stationary solution; see Figs 3A and 
4A for the corresponding numerical assessments. The specific form 
of the dynamics (3) and (4) already suggested that Ty = 0 is just a 
special case of ikl = T # 0. We therefore expected the existence of 
the solutions above not to be affected by introducing homoge- 
neous, finite delays %u = % # 0, even if this appeared somewhat 
counterintuitive. In fact, numerics confirmed this expectation as 
displayed in Figs 3B, 3C, 4B and 4C. 

The introduction of distributed delays xu instead of a single % 
changed results profoundly. To exemplify this, we distributed 
delays Tjy by drawing them at random from a uniform distribution 
over a certain interval. Recall that according to the transformation 
from (1) to (3), a distribution of delays T/y generally implies an 
equivalent distribution of phase shifts Ay. If Ay differed for all 
populations k and /, the stationary solution "P s tat(^i, ■ • ■ >^jv) °f the 
continuity equation (4) required the presence of many distinct 
phase clusters. We could prove the existence of that set and, 




Figure 6. Phase distributions in the presence of inhomogeneous coupling C'^" . A: Strong coupling (K 2 = 0.6). B: Weak coupling (A': = 0.2). 
As predicted by the phase derivation (9), inhomogeneity in the coupling matrix resulted in similar behavior as a distribution in delays; compare 
Figs 5A and 5C with panels Fig. 6A and 6B respectively. For reasonably weak coupling we found a widening of the phase distribution equivalent to 
the case z^l ~W(0,0.05) with homogeneous coupling. Again, an increase in coupling strength resulted in a concentration of 
doi:10.1371/journal.pcbi.1003736.g006 
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Figure 7. Neural mass coupling scheme for anatomical 



coupling C k ' 



For C[7 h) 



the SC matrix determines the coupling 



between excitatory units only. The SC matrix is the neuroanatomical 
coupling matrix; here three coupled excitatory/inhibitory pairs are 
shown. The other connections between excitatory and inhibitory 
masses serve to establish oscillatory behavior, resulting in a coupling 
matrix that is relatively sparse compared to C*° m) and Cff®. 
doi:10.1371/journal.pcbi.1003736.g007 



although the generic solution appeared similar to the homoge- 
neous delay case xu = t, it did contain M = N centroid values 6k 
instead of the small number M = 2 or M = 4 shown above. We 
depict examples of phase distributions for several parameter 
settings in Fig 5. Interestingly, the heterogeneity in Ty, or 
equivalently in Aid, agreed with weakening the between-popula- 
tion coupling K2 in that both cause a profound widening of the 
phase distributions; compare Fig 5A with 5C. That is, for a 
network with homogeneous structural connectivity, it is not the 
presence of delays per se that hinders synchronization but rather 
the distribution of delays (or the lack of coupling strength). 

Inhomogeneous coupling 

According to (3) both distributed phase shifts (or delays) and 
heterogeneous coupling may in principle result in inhomogeneity 
of phase coupling. In other words, distributed delays Tti and 
structural heterogeneity may yield inhomogeneity in functional 
connectivity. We therefore expected a heterogeneous coupling 
matrix with homogeneous delays %ki = T to be accompanied by 
desynchronization equivalent to the case of homogeneous coupling 
and distributed delays. To verify this, we used the inhomogeneous 
coupling sketched in Fig IB, which can be given more formally as 



-,(inh) 



(inh) 



,(inh) 



Winh) 



^11 



40 50 




- =1 

. 7~ Al'J.**. 




0.15 
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r; ; j • ?r 




0.12 



■ -'— r-"-"-.' 

. l ■_. l : l LM U i 



20 30 
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20 30 40 50 



Figure 8. Anatomical connectivity matrix and resulting functional connectivities as a function of delay distribution. A: Sparse 
neuroanatomical coupling matrix (4p serving as structural connectivity (connections are given in yellow). B: Structure-function correlation 
KQx >Pkl)- ^ : Ove ra 1 1 synchronizability <p<./>. Between-pair coupling strength K2 appears on the horizontal axis; while blue, red, and black lines 
correspond to (t^;) = 0,0.05, and 0.1 respectively. Note that an increase in delay distribution width increased synchronizability due to increased (xki}, 
but the accompanying higher xu variance decreased structure-function correspondence as predicted by (3). D-F: Corresponding spatial 
synchronization matrices p kl for delay distributions ta/ = 0, 1/(0,0.1) and ta/ ~£/(0,0.2) respectively (^2 = 0.1). Color coding is displayed in the 
most right panel; only excitatory nodes are shown. 
doi:10.1371/journal.pcbi.1003736.g008 
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with sub-population connectivities 



r (mh) _ 



-u 



_ _ r (inh) 



and 4™ h) = 



where ~W abbreviates ~W 



-3K 2 

N 



-U 



5K 2 



denoting that the off- 



(sub)diagonal entries were randomly drawn from a uniform 

K 2 

distribution centered around + — — . 

N 

The numerical simulations depicted in Fig 6 confirmed our 
hypothesis. For C^ nh ' and %kl = 1 = 0.05 we observed a widening of 
the phase distribution similar to that shown in Figs 5B and 5C 



where coupling was established by C ( 



{hom) 



but with delay 



distributions Tjy~W(0,0.2) and Tjy ~W(0,0.05) respectively. In- 
creasing coupling strength reduced the width of the phase 
distribution comparable to the switch from Fig 5C to 5A. 
Functional connectivity thus seems to result from an interplay 
between structural connectivity Cki and delay structure %kU 
Therefore both should be taken into account when studying 
functional connectivity in neuronal networks. 

Anatomical coupling 

The coupling matrices Cu considered so far were admittedly 
quite academic. However, these seminal examples did provide 
important insights that | as we will show here | generalize to more 
complicated and biologically plausible cases. We performed 
simulations using the coupling scheme displayed in Fig 7. The 

matrix cS** had the same structure as cS° m ' but now the blocks 



(N). 



C. 



Ik, -The 



were given by C^J' = SC, Cjj' = 0, and C- 
acronym SC stands for 'structural connectivity' that here refers to 
a neuroanatomical connection matrix as can be derived using 
DTI/DSI imaging [32,33] and /jc, is the identity matrix with K\ 
along the diagonal. To be precise, we used a binary form of the 
Hagmann connection matrix; see [24] for specifics of pre- 
processing. 

In line with earlier studies we quantified functional connectivity 
in terms of phase uniformity or phase locking value of the pair-wise 
relative phases, i.e. Pki = |<e ! ^'W - ^*W))|. Using this synchroniza- 
tion measure, simulation results can be best summarized in the 
form of functional connectivity matrices constructed from p kl 
values for all available pairs. The effects of delay structure T/y on 
these functional connectivity matrices are depicted in Figs 8D-F 
with the underlying structural connectivity given in Fig 8A. The 
functional connectivity matrix appeared rather sensitive for 
parameter values, as increasing coupling strength from ^2=0.1, 
which was the value used in Fig 8D-F, up to ^2 = 0.12 resulted in 
a fully synchronized network as can be seen in Fig 8C. The sudden 
synchronization is reminiscent of the phase transition towards the 
fully synchronized state at the critical coupling strength K c in the 
Kuramoto network. 

In a nutshell, from (4) we could deduce the mechanism 
responsible for the general finding that structural and functional 
connectivity are positively correlated [5,6]; see also Fig 8B. Our 
results clearly show that delay distribution affects both the spatial 
distribution of functional connectivity (Figs 8B; 8D-F) and the 



overall level of synchronization in the network (Fig 8C). The 
increase of overall synchronization is caused by a decreased phase 
shift A/a by which the phase dynamics (3) converges towards the 
Kuramoto model, i.e. the delay induces a change in stability of the 
(partially) synchronized state. 

Discussion 

We investigated the effect of time delays in the coupling 
between neural mass dynamics, where we consider an oscillatory 
regime, established by creating pairs of excitatory /inhibitory 
neural masses. Although we employed a specific neural mass 
model, we do consider our results generic because the mappings 
Cki^Dkl an d T£/->A,y are largely independent of the generating 
nodal dynamics, presuming that the time scales in the system are 
sufficiently separated; cf. Methods. 

By using this oscillatory dynamics to describe activity in certain 
brain areas, our approach links directly to the ongoing dispute 
about changes of functional connectivity in resting state networks 
(RSNs). There is growing evidence from experimental research 
that spontaneous brain activity during rest is highly structured into 
characteristic RSN patterns [1,3,34,35]. These activity patterns 
seem not to be the result of structural connectivity alone [5,36], 
but to reflect a non-trivial interplay between the neuroanatomical 
structure and dynamics [37]. The distribution of time delays 
involved in this dynamics may have an important role in shaping 
patterns of activity per se and neuronal synchronization in 
particular [15-17,38]. 

Key to our analysis was the reduction of a neural mass network 
to a system of phase oscillators summarized in (3). Several previous 
studies struggled with computational complexity when trying to 
unravel effects of delays vis-a-vis coupling on network dynamics 
[15-17,38]. By contrast, our analytic reduction 'readily' allowed 
for disentangling the contributions of both structural connectivity 
Cm and delays Xkl to the phase dynamics (3). Delays Xkl entered 
the phase dynamics as phase shifts Ay, given a proper time scale 
separation of oscillatory and phase dynamics. Furthermore, we 
found that heterogeneity in delays yields effects equivalent to those 
of heterogeneity in structural connectivitiy. That is, connectivity 
and delay effects cannot be easily distinguished when solely 
looking at functional connectivity patterns. 

The decrease of Q. as a function of delay, as depicted in Fig 2, 
agreed with the analytical findings of [13] as well as with our 
small-delay approximations oudined in the Methods. However, 
when further comparing the current results with the literature, one 
has to realize that some fundamental differences exist between 
general phase oscillator networks and our dynamics (3). One of 
those differences is the finite dimensionality of the system (3). We 
assumed every excitatory /inhibitory neural mass pair to represent 
a single brain area, by which the dimension of the system under 
study may be fairly low. On the contrary, most analytical work on 
phase oscillator networks considered the limit A->co [14,39-41] 
rendering one-to-one comparisons all but trivial. This can already 
be appreciated by the rather dramatic finite-size effects in 
Kuramoto networks [42-45]. Moreover, our structural connec- 
tivity C§ is rather atypical due to its strong + -asymmetry. 
Usually, the connectivity structure in similar phase oscillator 
networks comprises either fully homogeneous coupling or the 
entries Cki are distributed according to some unimodal distribu- 
tion [42,46]. An exception is [47], who investigated repulsive 
coupling, which is similar to the excitatory/ inhibitory connections 
in C^ om \ That study reported the presence of two anti-phase 
clusters reminiscent of the separate £7X-groups observed here. 
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The time scale separation in (3) and the resulting simplification 
of delays Ty as phase shifts Ay also hinders direct comparison with 
studies on delays in the Kuramoto network. This may indeed 
explain our seemingly contradicting results. For example, we did 
not observe the emergence of multi-stability mediated by specific 
(•^■^-combinations as reported in [12,48]. In those studies, regions 
in the (z,K) phase plane were found in which synchronization was 
entirely absent. This is clearly not the case in our study. We did 
find synchronized solutions irrespective of delays. This is not 
trivial, because for Ty = 0 the phase dynamics equation (3) may be 
reduced to a cosine-variant of the traditional Kuramoto model, 
which is known for its inability to display synchronized behavior 
[49]. By exploiting the pairing of the £/X-groups and the 
+ -asymmetry in cS° m ' we could map our averaged neural mass 
network to a conventional, fully homogeneous Kuramoto network 
via a mere transformation of variables. Hence our system can 
display synchronized behavior even for vanishing T values. This is 
consistent with [50] who did not find any qualitative effects of a 
phase shift a on the stability of the Kuramoto network. 

Apart from choosing random values as entries of the Cy-matrix, 
inhomogeneous coupling might also stem from creating distinc- 
tively different sub-populations in the network. It can then be 
studied by modifying within- versus between-network interactions. 
Particularly interesting in this respect is the occurrence of 
clustering in the network. When delays are not incorporated one 
needs either structural inhomogeneity [5 1] or higher-order Fourier 
harmonics, in combination with an appropriate phase shift [30] , to 
achieve clustering. The phenomenon of clustering is important in 
the light of the study of RSNs, i.e. the strong spatiotemporal 
organization observed in brain activity during resting state 
conditions [1,3,34,35]. 

We briefly consider a simple, low-dimensional example of three 
isolated excitatory nodes. We define a cluster as a number of 
(excitatory) masses attaining the same centroid phase value 6k- By 
denoting the excitatory and inhibitory centroid values as 9n £ and 
dyij., respectively and assuming ZkJc+N/2 = 0> from the phase 
derivation (3) and its corresponding continuity equation (4) one 
finds constraining equations 

m + txfiKi sin (6 Hl - 9 H£ + ^ 

( Ri 

aj8 — (D n sin A 12 + D n sin A 13 ) 

Rj 

aj8 — (Z>2i sin A 2 i + D 23 sin A 23 ) 

— * K£ 
Rj 

acp — (Z>3i sin A 31 + Z> 32 sin A 32 ) 

for 0 < Qx < —. In fact, these forms already hint at the interference 

between coupling and delays and its effect on synchronization 
structure. First, all terms on the right-hand side must have equal 
magnitude requiring specific combinations of D^t and Ay. 
However, both D^t and Ay are constrained by biology: Dy by 
the neuroanatomical coupling SC as part of Cy; and Ay by the 
spatial structure of the brain, as delays are proportional to distance 
between masses k and / due to finite conduction velocity. Second, 
because the left-hand side does not vanish, Z)y must have some 
lower bound. If Cy«l, then sin Ay cannot compensate this and 
the equality cannot be satisfied. Because Cy determines Dy, there 
must be some minimal coupling strength between nodes for 
synchronization to emerge. This explains the positive correlation 



between structural and functional connectivity, see, e.g., [6,17,38]. 
It also shows the intricate interplay between structure and delays 
in establishing synchronization structure. 

Interestingly, Z)y may be regarded as the effective coupling 
matrix that is typically encountered in dynamic causal modeling 
approaches [52]. The fact that Z)y is directly determined by Cy 
also explains the finding that models using the structural 
connectivity as a prior do show more evidence than models using 
other priors [53]. That is, models that have structural connectivity 
as a starting point, perform better in terms of data explanation. 
The sparsity of Dy induced by Cy may yield coexisting 
synchronized and desynchronized groups within the network, 
which are often labelled chimera states in the study of phase 
oscillator systems. It has been found that they crucially depend on 
the combination of coupling strengths and phase shifts [45,54,55] 
(or delays [56]), confirming that there has to be a specific matching 
of coupling and delays for synchronization to occur. 

Against the background of the aforementioned Ty-dependence 
of functional connectivity and the functional-structural connectiv- 
ity correlation for a biological plausible network, we numerically 
investigated this by performing simulations of (1) with structural 
connectivity given by the anatomical connectivity matrix 
reported by Hagmann and co-workers [32]. Functional connec- 
tivity was quantified as pair-wise phase uniformity, i.e. the phase 
locking value. Our numerical assessment is summarized in Fig 8. 
It clearly revealed off-diagonal patches with synchronization 
between nodes that are not coupled (contradicting what has been 
sketched above). The topology of the Hagmann et al. network 
shares similarities with the Watts and Strogatz' small-world 
network [57], i.e. both have a relatively large clustering coefficient 
with a small average path length. This kind of topology is often 
believed to be generic in biological neural networks like our brain 
[5] and enhances synchronizability compared to random networks 
[57-59]. The presence of sparsely connected clusters establishes 
synchronization between nodes that are only indirectiy coupled via 
their clusters. This causes 'blurring' of the structural connectivity 
matrix: The functional connectivity matrix is less sparse than the 
structural one [60] . Although this 'blurring' is similar to the effects 
attributed to volume conduction [61], in this case it is solely due to 
network topology. 

Next to such clustering phenomena, we can make even more 
general predictions about the effect of delays in this network. 
Structural and functional connectivity are most prominendy 
correlated for homogeneous delays, since Ay = A yields an 
interaction term in (3) that is merely a scaled version of C^p. 
Hence, the resulting spatial synchronization distribution largely 
resembles Z)y and thus presuming that the overall coupling 

strength is not excessively large. This effect can be seen in Fig 8B, 
where we depicted the Pearson correlation coefficient between the 
lower triangular parts of Cgp and the functional connectivity 
matrix p kl . Increasing the width of the delay distribution results in 
a decrease in structure-function correspondence. The positive 
correlations are consistent with the finding that the pattern of 
resting state activity is spanned by the eigenmodes of the 
underlying connectivity matrix in [62]. This is not as trivial as it 
may seem because the node dynamics in this study were noise- 
driven fluctuations around a stable fixed point and therefore 
entirely different from the self-sustained oscillations considered in 
the current study. 

Widening the delay distribution also had another effect: It 
increased its average value and consequently the mean phase shift 

71 

<Ay> = <— — Qty) tended to vanish. Therefore, the interaction 
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term Dki sin (<j>i — <j) k + A«) became more similar to an 'ordinary' 
sine-term, which is known for its capacity to enhance synchroniz- 
ability [63]. We illustrate this effect in Fig 8C, where overall 
synchrony <p> is shown as a function of coupling strength Ki for 
different average delay-values <T/a>. A similar phenomenon has 
been reported for a system of coupled Hindmarsh-Rose neurons, 
where a stable synchronized region appears to exist despite the 
presence of a (constant) delay [64] . 

Throughout this study we assumed the amplitudes to be 
constant. The relation between the envelope dynamics of M/ 
EEG and fMRI-BOLD signals [3] suggests that considering the 
temporal change of the amplitude may be very important for 
unravelling the spatio-temporal structure of resting state brain 
activity. Given that we focused on phase synchronization 
together with the slow time scale on which the BOLD 
dynamics evolve (<0.1 Hz, [1]), we believe that the assump- 
tion of constant amplitude is justified here. Investigating this 
assumption in depth, however, is beyond the scope of the 
current study. 

We summarize that the dynamics of a system of coupled 
Freeman neural masses (1) can be captured by the averaged 
phase dynamics (3), in which the role of the structural 
connectivity C k i and delay distribution x k \ become explicit. By 
this, one can identify the relative contributions of structure and 
delay to phase synchronization, i.e. to the functional connec- 
tivity of the neural network. Heterogeneity in structural 
coupling and distributed delays have equivalent effects on the 
observed phase distributions. Overall, this supports the notion 
that structure and delay are both crucial determinants of 
network behavior and should therefore be taken into account in 
unison whenever modeling realistic neural networks [37]. Our 
examples on clustering detailed how the intricate interplay 
between coupling and delays determines the form and spatial 
distribution of clustering in these networks. Pinpointing the 
explicit contributions of x k \ and C k \ in the phase dynamics (3) 
enabled us to understand their roles in establishing synchroni- 
zation structure and why functional and structural connectivity 
are so closely correlated. This implies that the observed 
temporal changes in synchronization structure in resting state 
and task conditions can be modulated through either x k i or 
amplitudes R k . 

Methods 

To analyze the gross membrane voltage of a neural mass we first 
wrote its dynamics (1) as a two-dimensional system using the 
auxiliary variable Uk=Vk- For the sake of simplicity we further 
rewrote the convolution integral in (1) by means of 
[W(t) * S[V,{t-z k ,)}} = [W{t-z k i) * S[V,(t)]]. Then the dynam- 
ics (1) reads 

Vk=U k 

u k =-(<x k +P k )u k -<x k p k v k +<x k p k q k ... 

(0) 

N 

+ x k Pk 52c u [w{t- x u ) * S[ V,(t)]] . 
1=1 

In the following we discuss this dynamics in its oscillatory 
regime after transforming the system into polar coordinates to 
derive the corresponding phase dynamics. That transform, 
however, was not applied to the original state variables [V k ,Uk] 
but to the deviations [SVk,SUk] around the unstable fixed points 



[F< 0) ,[/< 0) ], i.e. we transformed [V k , U k ] -> \vf > + SV k , Uf ' + 

5U k ]. Furthermore, we expanded around in the 

sense of Taylor and obtained 

s[vf^ + 8V k {t)] = Y J ^[vf ) ]{6V k ) n ■ (7) 

With this expansion the dynamics of the deviations 5 V k ,SUk could 
be approximated as 

sv k =su k 

5U k n-(a k +p k )SU k -* k p k dV k ... 

W 

N 

+ ockpk J2c k ,\w(t- x u ) * S' \ Vf ] ] S K/(f)l . 
/=i 

A closer look at this linear system revealed that the fixed points 

were indeed unstable nodes, provided that a proper 

balance between excitatory and inhibitory masses was present. 
That is, the eigenvalues of the linear system (8) were complex- 
valued with positive real-parts — explicit expressions for the 
eigenvalues (as a function of delay) can be found below. 

Next, we transformed system (8) into polar coordinates by 
means of [SV k ,SU k ] ->[R k cos(Clt + <l> k ),-ClR k sm(Qt + <j) k )] so 
that the phases could be cast into the generic dynamics 

h=-n+-^(zU k SV k -5V k dU k ) . (9) 

Importandy, in the current case we could assume that the phase 
dynamics (9) in general contains two (or more) distinct time scales: 
the rate of change given by the oscillation defined via the 
frequency £1 and the time scale(s) of the amplitudes Rk and the 
phases <j> k ; the latter are much slower than In/Q, and can hence be 
separated. More formally we used 

\R k /R k \«Q. and |^/^|«Q . 

In first order approximation we could thus consider R k and <f> k to 
be constant during one period of oscillation, T = 2n/Q. This 
approach is referred to as a combination of a rotating wave 
approximation and a slowly varying amplitude approximation [65] . 
It enabled us to average the dynamics over the interval [0,2"); see 
also [23]. As will be shown in the following, this averaging 
procedure decoupled amplitude and phase dynamics, which 
ultimately resulted in the dynamics (3). Below we will provide a 
more formal discussion regarding these approximations. 

Averaging — towards the phase oscillator model 

To average the dynamics (9) we defined the averaging operator 
as (f{s)}= \ lf(s)ds. We first substituted (8) into (9) and used 

<5V k 5U k } = 0, <<$^>=^l,and (SUf c }=^ . 

The convolution integral on the right-hand side of the second 
equation of (8) required more attention. Recall the definition of W 
in (2) and the definition of the convolution operator, with which 
one can write 
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[W{t-x kl ) * SV,(t)]= W(t-x kl -s)8V,(s)ds 
= | -co ye-t'-W-^Ri cos (£h+<j>,)ds 



= R,ye->i'- T ki) 
yR, 



z kl „ys 



e ys cos(Q.s + </>/) ds 



y 2 +n 2 



{ycos(Q(f-T W )+^)+£2sin(Q(f-T H )+M} .(10) 



When multiplied by SV k (t), this yielded the two averaged 
trigonometric expressions 

<^cos(£2f + ^)^^cos(Q(f-T H )+^)> 



1 y 2 RiRk 

2y 2 + fl 2 



{cos (fiTw) cos(0; - <t> k ) + sin (fliw) sm(<j> l - <f) k ) } 



and 



Q.yR, 

y+o 2 



</Jtcos(£2f+^)^- 7 ^sin(£2(f-T«) + ^)> 



1 flyj?/^ 

2 y 2 + fl 2 



{cos (Qtjh) sin(f - <jy - sin (Qt«) cos(<^, -<j> k )} 



After substituting this in (9) we obtained the phase dynamics 



h = a k - a k fi k J2 A kl sin (</>, - <f> k ) - 
1=1 

N 



(11) 



i=i 



where we defined the constants A k \ and B k i as 



y 2 sin(nTA-/) + nycos(nT /rf ) R, r (0) 



2Q(y 2 + n 2 ) 



R k 



B k i = 



y 1 cos (Qt«) - £2y sin (Qt«) 



2Q(y 2 + n 2 ) 



[>f]c„ 
[>f]c„. 



By this procedure we omitted all fast oscillating terms as they 
averaged out (cf. rotating wave approximation). 

As mentioned in the Results section, we focused on the case in 
which our convolution kernel W did not contain any memory. 
That is, we considered the limit y— >oo. In this limit only the first 
terms in the numerators of A k \ and B k i remained non-zero and we 
could cast (11) in the form (3) using 

D kl = ^Al l + B 2 kl = >] C« and A kl = * -£lx kl . 



Note that in this form, the delays x k i only appeared in the phase 
shift A k i. Last but not least we simplified expression (3) by 
exploiting the homogeneity of CjJ om ' to explicitiy formulate the 



D k i matrix multiplication. In particular for equal delays, i.e. for 
x k i = x, this led to a greatly simplified form of (4) that we 
summarized in the Results section and will be discussed in more 
detailed below. 

Fixed points and amplitudes 

The neural masses do not oscillate around the origin but around 
the fixed points ^Pjt"\Cj^J , which have a direct influence on the 

coupling terms D k [ through the term 5" ^ I 7 / 0 ' J . Delay values x k i do 

not influence the positions of the fixed points [Pj^.t^J because 

by definition [j^C/fl =[0,0]. Hence, V k {t-x kl )= V k (t) holds, 



presuming the fixed points exist. Therefore we were free to choose 
x k / = Q, such that under the limit y->co the coupling term in (6) 
reduced to 



f r i r i 

a k fa £ C U 1 X W(t-x kl -s)S\ Vj 0) ds^a k fa £ C k ,S\ vj 0) 

1=1 J 1=1 



After inserting the form of C^° m * we explicitly found 



yf> = 

l<tk 1 1 l^k + N/2 11 1 1 



._«/•> L J L J L J 



el 



I?k 
lex 



which implied Vf ] =vf\ Vief and vf> = Vf\ V/ceX. 
That is, in the case of the homogeneous coupling the fixed points 
of the excitatory masses are equal and the same holds for the 
inhibitory masses. 

The coupling D k i also depended on the amplitudes R k of the 
neural masses — see also [24] . Accounting for the high degree of 
homogeneity in the system, we assumed the amplitudes to be equal 
for equal types of neural masses, i.e. R k = R£, k e £ and R k = Rj, 
keT. Furthermore we randomized the parameters at = a + e, 
P k = P + e by introducing £ as a mean-centered random variable. 
Whenever appropriate we chose £ sufficiendy small to restrict 
discussion to the mean values a=\/N ~}2 k= j 0l k and 
/?=l/^£Li Ac- 
Homogeneous coupling — existence of solutions 

Since the phase oscillator system (3) can be cast in Kuramoto 
form, fully synchronized solutions may be stable despite the 
presence of equal delays x k \ = x. But how about solutions other 
than the fully synchronized ones? In what follows we discuss 
existence and form of partially synchronized solutions of (3) 
for general delays x k /. We concentrated on homogeneous 
coupling and varied the distribution of x k \. In the homogeneous 
case we found the dynamics of the &-th node's phase distribution 
V k to be 



j t v ^,t) = --^- |Q i ,-ti'Y.M+> <i':,+\'.i)r l {,i>,.t)M ! 



- Hz r E «n U, -h + A H )V, U„t) \V k {4> k ,t) 



(12) 



where the subscript ° = £ when k e £ and ° = T when keT. Note 
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that equation (12) may differ for every k. Homogeneity of Cj^ om ' 
enable us to express Dki explicitly and to define the following 
constants 



Tei = a/? 
T x£ = aji 
Tjj = apK 2 



Rx 
Re 

Rs 
Rx 



vf\+K 2 - 



N-2 



Vf\+K 2 - 



N-2 



N 



(0) 



-S \ V; 



(") 



(13) 



Sufficient for the existence of a stationary solution is the case in 
which the drift coefficient in the dynamics of the probability 
density vanishes, here the bracketed term on the right-hand side of 
(12). From the dynamics (12) it readily follows that the phase 
distribution obeys the form (5), i.e. V = V static \, ■ ■ ■ An)> contain- 
ing N different centroid phase values 8k- 

Equal delays. Ty = T. As said, for homogeneous coupling the 
sums in (12) could be evaluated in the form of the constants in (13). 
Furthermore in the case of equal delays Ty = * (i.e. Ay = A) we 
only had to account for two distinct populations each with just a 
single type of density dynamics, by which the system (12) of N 
equations could be reduced to merely two distinct ones: 



fs L 



co — Tss 
■In 



r-<j> e + A)VeU e> t)d(l>g 



dt 



co-T xe sinU £ -fe + A )VeUe,t) 



JO 



(14) 



To show the existence of 'P s tat,2(^£^x) consisting of two phase 
clusters, one for the excitatory units and one for the inhibitory 
units, i.e. 

7W0 £ A) =\&(hSe)+\ Kti - Ox) 

we followed a constructive approach and substituted Vsts>t,2(<l>£,<!>x) 
into (14). Vanishing of the drift coefficient required 

{Tee + Txx) sin A = T £Z sin(A + xj/) + Y JZ sin(A -i/,). (15) 

when abbreviating \j/ = 8 £ — 8j. Using \Tx£ + Tsx\>Wes + ^xx\, 
which followed from (13) irrespective of the value of K 2 , we could 
conclude that a solution \f/ satisfying (15) exists. Note that only i/f 
appeared in (15) and not the centroid values 0„, which allows for the 
mapping to the conventional Kuramoto model as will be discussed 
below. 

We could readily generalize this line of reasoning to an arbitrary 
number of clusters after defining the stationary phase probability 
density containing M clusters as 



, M/2 . m 

p sta t,M(0 1 ,...i M )=T7E ,5 ^- e ^)+M E Wi-h*) 

k=\ M k = M/2+l 



We defined generalized phase differences i/^yu = Q a ,k — 6b,I 
where a,b e {£,1} and k,l e {1, . . . ,M}, by which we obtained 
the set of constraining equations as 



E sin (^££.« + AA/ ) 



+r IX 
=r £ j 



Isl 

E sin (^.« +A/rf ) 

JeX 

E sin (^.« +A/rf ) 



(16) 



leE 



In contrast to (15) we here considered general delays tkh 
but note that the form of (16) did not change for T&/ = T. As 
the values of the constants T , still fulfilled the inequality 
|Fj £ + r £ i| > |T ££ + Till we again could conclude that this system 
can be solved. Thus, even for equal delays %kl = t (or Ay = A) the 
stationary phase distribution may contain any number of clusters 
M <N. Of course, this existence proof does not imply that all 
these solutions will be found in reality since we have not yet 
addressed their stability; cf. Figs 4A and 3A. Note that apart from 
the fully synchronized state, we here restricted our stability analysis 
to numerical evaluations. 

The distributed delay case (16) does not require Ay = A. Hence, 
for general Ay the stationary distribution V s tnt,M will contain 
M = N centroid values 6t with k 6 {1, . . . M}. 

Stability for Ty = T — a fling with the Kuramoto 

network. Before proving the existence of distinct solutions in 

the case of arbitrary delays Ty, we first briefly discuss the case of 

constant delay because it readily links the model (3) to the 

Kuramoto network of coupled phase oscillators [63]. For 

Ty = T = 0, the interaction term in (3) reduces to a mere cosine 
71 

term because of Ay = — . We note that in Kuramoto's traditional 

form this cosine term leads to a system that is unable to show 
synchronized solutions [49] . Our numerical simulations, however, 
revealed synchronized solutions; see, e.g., Figs 3 and 4. In what 
follows we show that due to the fact that the centroid phase 
values Ok are irrelevant one can distill the traditional Kuramoto 
model. 

It is the coupling matrix C^ 0 ™' that made the difference with 
Kuramoto's network as we assumed a (balanced) combination of 
excitatory and inhibitory units. In the given form illustrated in 

Fig 1A, the matrix C^° m * is strictly speaking not homogeneous but 
contains an asymmetry. This, however, did allow for mapping our 
model to the Kuramoto network by adapting an approach of 
Frank and co-workers [31]. In essence, we applied a change in 
variables yielding a fully homogeneous coupling matrix. For the 
sake of legibility we set R £ = Rj, S'[vj 0) ] = S' [ V^] and K x = 2K 2 . 

Then Dy simplified to Z)y = + ^'°'] wnere the ' + '-sign 

discriminates rows that correspond to the excitatory /inhibitory 
neural masses £/T, respectively. By this simplification the key 
feature of the system (9), the asymmetry in the coupling matrix, 
remained untouched. 

We further transformed the system into a rotating frame at 

a/?-n 2 

frequency co = — — — by defining [65,66] 
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<Pk : 



<j> k — wt, ke£ 
j> k — 7i — cot k el 



(17) 



After abbreviating the constants ^A^S" [F (0 '] =jc, the phase 
dynamics (9) became 



In line with Kuramoto's analysis we introduced a mean (or 
cluster) phase 9 and amplitude p in terms of 
,a 1 x-^JV 



pe 



N 



e ''^' +A ' that represent well-known order param- 



eter(s). Finally, substituting 

^ i sin (^ ~ ( Pk + A )=P sin ( 6 ~ <P) yielded 



ip= —Kp sin(i9 — (p). 



Here we dropped the subscript k because the equation applies for 
all nodes. This dynamics agrees entirely with the Kuramoto model, 
apart from the fact that we here considered cci c = oc and fa = by 
which all oscillators' natural frequencies w k = w agree. Due to this 
correspondence to the Kuramoto model, we could conclude that the 
system with equal delays Tkl = T and 'homogeneous' coupling matrix 
CjJ om ' can generate synchronized states provided the overall 

coupling strength — K2S' [F (0 *] is properly chosen. 

Numerical assessments 

Both distributed delays Tkl and heterogeneous coupling called 
for numerical assessments, particularly when it came to the 
stability of solutions. We performed numerical simulations of the 
coupled neural masses (1) using a conventional Euler-forward 
integration scheme with time step At= 10~ 3 .s over a time span of 
10 s. We verified the appropriateness of this simple implementa- 
tion against a more elaborate predictor/corrector integrator [67], 
which revealed little to no difference but demanded far more 
numerical resources. The simulated network consisted of 500 
nodes (250 £/X pairs) with a k and fa being randomly drawn from 
uniform distributions (0C4 6 [25,75] and fa e [175,225]) to mimic 
distributed natural frequencies per £ /X-pair of nodes. Although 
randomly drawn, these sets were fixed across simulations trials. 
Initial conditions where chosen randomly and did vary between 
trials; V k (t = Q) e [- 100, 100], A: e £, V k (t = 0) e [-W,60],k e 1 
and similarly for Vk(t = 0). The external input qk was set to 
qk =20, ke£, q/ c = 0, keX. Coupling between masses was 
achieved by using the sigmoidal activation function S[ V] that was 
given as 



S[V\=y 



1 



\ +e -{y-y&)i° 



whereas the constant y is just a scaling factor [68]. In the 
simulations we used the following values for the activation function 
parameters: y = 250, K t h = 15 and <r=l. 

In order to compare the numerics with our analytical results, we 
determined the phase values <j) k from the simulated potentials Vk{i) 
by means of the Hilbert phase. To this end, we first determined 
oscillation frequency £1 as the lowest frequency with a coinciding 
peak in the power spectra for all nodes. Voltage traces Vk(t) were 
band-pass filtered using a 1 -st order bi-directional Butterworth filter 
in the frequency band [0.8fl,1.2Q]. For each sample in the interval 
t e [7,8], phase values (j> k (t) were then calculated as the angle of the 
analytical signal. By restricting analysis to that interval we avoided 
transient behavior as well as possible filter artifacts. The so- 
determined (j) k contained the frequency component Qt, which we 
first subtracted to obtain <p k . Then, we opted to compute phase 
distributions over the time interval T and over successive trials. The 
(circular) mean phase, however, differed from trial to trial because of 
the randomly chosen initial conditions (see above). Hence for every 
trial we shifted phases by the mean phase of the excitatory 
population prior to concatenating trials. By this, the phase 
distributions of the excitatory phases always became centered 
around zero and the inhibitory phases were considered relative 
values. The mean phase per trial was given as 

<D £ = arctan' ^ ^ sin q> k {t)/ ^ ^ cos Vk (t) 

\ktEteT kcEteT 

where arctan' denotes the quadrant-corrected inverse tangents. 

The distributions displayed in the figures are the phase 
distributions obtained from 100 simulated trials with different initial 
conditions. As said, parameter values were identical across trials. 
For the simulations involving C?p connectivity, we used a binary 
form of the 66 areas parcelated Hagmann et al. matrix [24,32] as 
C^P block; see Fig 7 and 8A for the coupling scheme and Cgg\ 
respectively. Functional connectivity was quantified as phase 
coherence given by p k / = |<^e'('W')~'M')))>|. We performed one 
hundred simulation runs of 10 s for each Tyy distribution with 
different initial conditions and a.k,fa and averaged the p k i matrices 
over these runs. This was done to avoid high p kl values due to 
common oscillation frequency alone. For each run, data in the 
interval t e [4,8] was used to determine p k {. The overall coupling 
strengths were set to K\ = 1 and K2 =0.1 for the p kl matrices 

displayed in Fig 8. Structure-function correspondence r^C'^p ,p k i^j 

was quantified as the Pearson correlation coefficient between the 
lower triangular parts of both matrices to avoid spurious correlation 
values due to common terms along the diagonal. 

Eigenvalues of the linear system 

To estimate the oscillatory regime of the system of coupled 
Freeman neural masses (1) we considered the linearized dynamics 
(8) for y — ►co which in general reads 

SV k =5U k 

SU k « - (a* + fa)SU k - XkfaS V k 

N r 1 

+ a k faY J C kl S l [vf ) \6V,(t-T kl ) . 



The fraction in this equation may be interpreted as the 
cumulative distribution function of the normal distribution 
J\f(V— K t h,ff 2 ) of the firing thresholds V t ^ across the population, 



Assuming all delays T k i to be (very) small we expanded 
5Vi(t — z k i) in the sense of Taylor and approximated up to the 
linear order in Tkl- 
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SV,(t-T U )xSV,(i)- th j ( i V,(t)} =dV,(t)-T kl 5 U,(t). 

For the sake of legibility we here considered a single isolated £/I- 
pair with Cge = Cxi = 0 an d Cxs = — Cfz = K\ . We also assumed 
equal delays, i.e. x k i = x. Then, we found the resulting linear 
dynamics as u = AM where we abbreviated u = 
(dVc,dU £ ,dV x ,dU x ) T and 



A = 
( 



0 1 0 0 \ 

— ajS -(=< + /*) -a/?^i5'[l4 0) ] ^jSK i S'\vf ) 

0 0 0 1 

oc/?Jfiy[F* 0) ] -t«^iS"[)^ 0) ] — a/S -(a+j8) y 



The matrix ^4 came with eigenvalues A,- 



Q(t = 0) = 



1 

2a/2 



(a - ft 4 + 1 6ot 2 fi 2 Kj S' [ K< 0) ] 5' [ Ff] - (a - ^) 2 



Further, for the real-part to be positive the coupling constant K\ 
had to be sufficiently large and the intrinsic damping a and/ or f} 
sufficiently small but finite, because of 



»{A!, 2 } T=0 >0 oK,> 



where we note that 3R{A3 4} < 0. 



zpS>[vf}s'[vr} 



y ' + 2i^[(<x + /J)T-2A: 1 ] + [(a- p) 2 -Kf^} 



where we abbreviated \jj = oifi^J S' [V^'J S' |Vj 0) ] . These eige 
values had the following real and imaginary parts 



2 1 



\J [(« + Pf -Kj^ 2 T 2 ] 2 + # 2 \{a + P)x- 2Ki} 2 +[( a + P) 1 -Klfx^ 



^\j\J[(v. + P?- Kf ,/, 2 t 2 ] 2 + # 2 [(a + /?) t - 2JC, ] 2 - [(a + j^-XftrV] 



A necessary condition for the existence of a stable limit cycle, 
and hence for the system (1) to display oscillatory behavior, is that 

the fixed point {^Vf\uf\V^\Vj^ is unstable. This means that 

for at least one of the conjugate pairs Xj, Si(Ay) > 0 must hold, 
which was indeed the case irrespective of x. The corresponding 
Sr(A/) then provided a rough estimate for the frequency fl as a 
function oft, as shown in Fig 2 (solid line). In the particular case of 
t = 0 we found 



Distinct time scales 

For the separation of time scales underlying all the major 
approximations in the current study we considered the case in 
which two distinct time scales are present in the system of coupled 
neural masses: the oscillation described by the (mean) frequency Q 
and its corresponding period 2n/Q — from here-on referred to as 
fast time scale — as well as a slower time scale, on which the 
dynamics of <j) k evolve. In what follows we will verify the 
expression for <j> k in (9) and show how the separation of time 
scales enabled us to determine the role of x k i 10 the convolution 
W(t — x k l) * S[V[(t)]. As is conventional in multiple-scaling 
approaches, we set the time t as 'fast' time and the 'slow' time 
as r\ = et with 0<e«l. For the sake of legibility we here adopted 
the dot-notation for temporal derivatives d/dt and further 

d 

abbreviated partial derivatives as =Oq. 

We denoted the deviation of the voltage from its fixed-point as 
5 V k (t) = Rk cos (Qt + <f> k ) where we assumed that <j> evolved on the 
slow time rj, i.e. ^> k = ^ k {f])- Note that an equivalent approach can 
be adopted for the amplitude dynamics, i.e. Rk = Rkil), which is 
referred to as the slowly varying amplitude approximation. 
As we were primarily interested in the (j>k dynamics, we 
regarded amplitude Rk as constant on both time scales; see 
Discussion section. By this we could readily apply the chain rule 
and obtained 



dV k = d,dV k + £d v dV k .= -(Q+ed^ k )R k sia(Qt+^ k ) = dU k 

where the last equality follows from (8). For the derivative SU k we 
found 



»{4/} t _ 0 =-^±^\l^(x-P) 4 +16Kfy 2 + (a-p) 2 SU k =-(Q + edJ k ) 2 R k cos(Qt + ^ k )-e 2 8 2 ^ k R k sm(nt+^ k ) 



0 ± 2y/2 



1 ^(x-tf + l^-ia-P) 2 



revealing that for K\ #0 the imaginary parts did not vanish, i.e. 
the £/2T-unit always displayed oscillations around the fixed point 
[F^ 0 ', Fj 0 ^] because the sigmoid's derivative is positive definite: 



Next we considered the expression (9) for <f> k , where we 
emphasize that § k could here be identified as d n $ k , since <f> k 
evolves only on the slow time scale. To anticipate: (9) discarded all 
terms evolving on the fast time scale, i.e. all 0(1) terms in favor of 
0(e) terms and higher. To show this, recall that the right-hand 
side of (9) read 
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Cl+--^(dU k dV k -8V k dU k ) 

UK 



= -Q + 



{a 2 + 2rf»„fc + e 2 d 2 n h)Rl ~ ^ e 2 ^ k R\ sin (2Qt + 2<j, k ) 



2ed v <j> k + £ - d 2 ^ k (l - X - sin (2Qt + 2<j> k ) 



In words, only expressions evolving on the slow time scale were 
retained, i.e. only O(e) and 0(e 2 )-order terms. When focusing on 
the slow time sale 0(e) and discarding the even slower time scale 
0(e 2 ), we could conclude that, up to a constant, (j) k is given by (9); 
note that we here applied the so-called two-timing method; see, 
e.g., [23,69,70]. 

As said, we used the constancy of <j> k (on the fast time scale t) to 
evaluate the convolution term W(t — x k i) * S^vj 0 ^ . We exploited 

the description in two time scales to justify the transformation of 
the delay x k \ into phase shifts Qx k i. We explicidy evaluated the 
integral e' /s cos (Qs + ^^ds to show that (10) is its 0(1) 

result. For the sake of readability we dropped the explicit time 
dependence of <f>i whenever possible. 

First, by integrating by parts twice we obtained 



e^cos (Qs + ^^ds 



= -e y ''cos (Qs + 0,) 

Y 



1 f'~ T « 

j (Cl+sd^,) sin {Q.s + ^,)e ys ds 

J — CO 

- e' { '-'k,) cos (Q(? - t«) + U<t ~ *»))) 



+ e ys {Q. + edj,) sin (as + <j>,) 

yZ 



-*ki 



1 t'-^kl , 

-2 (Qs+e5^,) 2 cos(ay+^,)e y '' 

J J — CO 



+ e 2 d 2 J, sin (Qs + <f>, )e ys ds . 



Then, when discarding 0(e 2 ) terms, we found that the integral 
with which we started appeared again on the right-hand side. This 
allowed us to write 



~ I [e^'-w) cos (fl(r-TH) + -t«))) + 
- (Q+ ed^,) sin (Cl(t - x kl ) + <j>,(e(t - t w )))] , 



which resulted in 



e' s cos (Qs + ^yds 

- CO 

e*-'ki) [ cos (tl(t - t„) + " *«))) 



y 2 + n 2 +2e^ / 

+ (Q + ed^i) sin (Q(/ - x kl ) + $,(e(t - x k ,)))] . 

With this form we could finally express the convolution terms as 



W(t-x k ,)*dV,(t)- 



yR, 



[y cos (; 



0.{t-x k ,) 



'y 2 + Q 2 + 2ed^ l 

+ <t>Mt- x k ,)j) + (O + ed n h) sin (d(t -x k ,) + <j>,(e(t - x k ,)))] . 

For e « 1 time scales are sufficiendy separated to ignore all 0(e) 
terms by which we arrived at (10). 

The derivation of (10) required the intuitive assumption 
<l>l(t — x k i)x<l>i(t), which might be motivated by a (relatively) 
small-delay approximation. Consider the Taylor expansion of 
<j>l(t — x k i) around t, which reads 



fa(t-v u )K4,(t)-—[4 l (t)}x kl = 4,(t)-ed^ l (tyz kl . 
If e«l and x k i is of the order 0(1), i.e. of the same order of 

271 

magnitude as the oscillatory period T= — , then one may 

conclude that the approximation (j>i(t — x k i)x(j)i(t) is valid. Note 
that this consistent with [71]. In particular, when x k \ is of order 
0(1 /e), i.e. of the same order of magnitude as the slow time scale, 
the approximation (j>i(t)~<j>i(t — x k i) fails. 
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y 2 + n 2 + 2gg^ / 



e y,! cos (Qs+fi^ds 
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